#ifndef DDAS_MAIN_H
#define DDAS_MAIN_H

void updatePercentageDone( verylong current_result_count, verylong total_result_count );

void ddas_worker(gsl_matrix* independent_vars, gsl_vector* coeffs, gsl_matrix* covariance, double *chi_squared);


/*
* @author Max Rupplin
* @date 5-15-08
* 
* 'worker' is an implementation of BOINC's 'worker interface' and is used as a function pointer for creating the 'worker' thread via 'boinc_init_graphics(worker)'
*/
void worker()
{
	//holds covariance results for a single result
	gsl_matrix* covariance;

	//holds the 'on' columns for a given iteration ('101' -> 1st, 3rd indices are 'on')
	gsl_matrix* required_independent_variable_columns;

	//holds the regression coefficients for a single result
	gsl_vector* regression_coefficients;

	//how many results have been processed so far, relative to starting_value (if starting_value is 5000 and iteration=5001, current_result_count = 2) 
	verylong current_result_count=0;	

	//simple for-loop iterator
	verylong i=0;

	//the SSE (sum square error) for a single result
	double chi_squared;

	//set % of work done to 0
	boinc_fraction_done(0.0);

	//read in data for work scoping - this populates or informs functions to populate, the global variables found in ddas_globals.C
	ddas_parse_scoping_file();

	//populate both independent & dependent data into a single gsl_matrix, 'master'
	ddas_parse_data_file();


	//using 'master' populate 'independent' matrix
	ddas_populate_independent_matrix();

	//using 'master' populate 'dependent' vector
	ddas_populate_dependent_vector();

	//starting_value and result_set are potentially defined by state file (if previous state exists)
	if( !ddas_parse_state_file( &i ) )
	{
		//since the state file was empty, i = starting_value (just start at the beginning)
		zcopy( starting_value, &i );
	}
	else
	{
		//current_result_count = i - starting_value
		zsub( i, starting_value, &current_result_count );
	}

	//a basic sanity check
	if( zcompare( i, starting_value)==-1 || zcompare( i, ending_value )==1 )
	{
		cerr<<"Error: ddas_main: i < starting_value || i > ending_value"<<endl;

		return;
	}

	//for all binary combinations between 'starting_value' and 'ending_value' perform gsl_multifit_linear and store top 'results_to_return' in 'result_set'
	for(; zcompare(i, ending_value)!=1l; zsadd(i,1l,&i))
	{	
		//update how many results have been considered
		zsadd(current_result_count,1l,&current_result_count);

		//how many '1' columns in binary representation of current iteration 'i'
		columns_used=ddas_get_one_columns_count(i);

		//get '1' columns from binary representation of current iteration 'i'
		required_independent_variable_columns = ddas_get_appropriate_independent_columns(i);	

		//allocate memory - calloc zeroes values
		regression_coefficients = gsl_vector_calloc(columns_used);

		//allocate memory - calloc zeroes values
		covariance = gsl_matrix_calloc(columns_used, columns_used);

		//call gsl_multifit to solve
		ddas_worker(required_independent_variable_columns, regression_coefficients, covariance, &chi_squared); 

		//consider storing new result
		if( ddas_get_max_SSE(result_set) > chi_squared)
		{
			#ifdef DEBUG_MODE
			cerr<<"A new best SSE was found."<<endl;
			cerr<<"The new value is: "<<chi_squared<<endl;
			cerr<<"The old value is: "<<ddas_get_max_SSE(result_set)<<endl;
			ddas_print_statistical_data(regression_coefficients, covariance, &chi_squared, &i);
			#endif

			//create new result
			DDAS_CLIENT_RESULT* tmp = ddas_client_result_alloc(i, regression_coefficients, chi_squared);	

			//replace worst with new (value replacement)
			ddas_replace_worst_result(tmp);

			//free memory
			delete tmp;			
		}

		//update boinc_fraction_done
		updatePercentageDone( current_result_count, total_result_count );

		//free memory
		gsl_matrix_free(required_independent_variable_columns);
		gsl_matrix_free(covariance);
		gsl_vector_free(regression_coefficients);

		//every 10,000th iteration check if it's time to checkpoint, if so do so
		if( boinc_time_to_checkpoint() )
		{
			//save state
			ddas_write_state_to_file(i);

			//notify BOINC that the checkpointing is complete
			boinc_checkpoint_completed();
		}
	}

	//rewrite state output
	//rewrite final results

	#ifdef DEBUG_MODE
	ddas_print_client_result_set(result_set);
	#endif

	//serialize results
	ddas_write_results_to_file(result_set);

	//free global variables
	free_globals();

	//free verylongs
	zfree( &current_result_count );
	zfree( &i );

	//tell BOINC computation is complete
	boinc_finish(0);
}


/*
* @author Max Rupplin
* @date 5-15-08
* @param independent_cols The independent data columns
* @param coeffs The gsl_vector used to store coefficient values
* @param covariance The gsl_matrix used to store the variance-covariance values
* @param chi_squared The SSE of the resulting computation
*
* 'ddas_worker' calls gsl_multifit_linear on a set of independent and dependent variables and stores the result in coeffs, covariance and chi_squard
*/
void ddas_worker(gsl_matrix* independent_cols, gsl_vector *coeffs, gsl_matrix *covariance, double *chi_squared)
{
	#ifdef DEBUG_MODE
	cerr<<"ddas_worker called."<<endl;
	//cerr<<"width of workspace: "<<columns_used<<endl;
	//cerr<<"height of workspace: "<<row_count<<endl;
	#endif

	//memory allocation
	gsl_multifit_linear_workspace* work = gsl_multifit_linear_alloc(row_count, columns_used);

	#ifdef DEBUG_MODE
	cerr<<"covariance vector before:"<<endl;
	ddas_print_gsl_matrix(covariance);
	cerr<<"independent_cols before:"<<endl;
	ddas_print_gsl_matrix(independent_cols);
	#endif

	//do linear fitting
	gsl_multifit_linear(independent_cols, dependent, coeffs, covariance, chi_squared, work);

	//free memory
	gsl_multifit_linear_free(work);

	#ifdef DEBUG_MODE
	cerr<<"covariance vector before:"<<endl;
	ddas_print_gsl_matrix(covariance);
	cerr<<"independent_cols after:"<<endl;
	ddas_print_gsl_matrix(independent_cols);
	cerr<<"ddas_worker returning."<<endl;
	#endif
}



/*
* @author Max Rupplin
* @date 5-15-08
* @param current_result_count The number of results already calculated at this point in time.
* @param total_result_count The number of results that will be calculated when this result is done. (ending_value - starting_value)
* updatePercentageDone is a helper function used to handle the calculation the input to BOINC's boinc_fraction_done(double).
* It can also handle cases where its inputs are numbers (much) larger than the largest doubles, in which case it uses some radix manipulation with LIP.
*/
void updatePercentageDone( verylong current_result_count, verylong total_result_count )
{
		//maximum long value (2^63)
		verylong max_long = 0;
		zsread("9223372036854775807", &max_long);	

		//do a rough approximation (cast to double) if the numbers aren't too big
		if( zcompare(total_result_count, max_long)==-1 )
		{
			//cout<<"crc: "; zwrite( current_result_count ); cout<<"trc: "; zwriteln( total_result_count );
			boinc_fraction_done( ((double)ztoint(current_result_count)) / ((double)ztoint(total_result_count)) );
		}
		//we are dealing with very large numbers
		else
		{
			//we'd like at least 4 digits of precision (75.01)
			int precision = 4;

			char crc_c[128];
			char trc_c[128];
			
			zswrite( crc_c, current_result_count );
			zswrite( trc_c, total_result_count );
			
			unsigned int crc_len = (unsigned int)strlen(crc_c);
			unsigned int trc_len = (unsigned int)strlen(trc_c);
			unsigned int ul_len  = (unsigned int)strlen("9223372036854775807");

			int diff = trc_len - crc_len;

			//need to think this thru - can diff-precision <= 0?
			if( crc_len <= 1 || ((unsigned int)(diff-precision)) > ul_len )
			{
				boinc_fraction_done( 0.0 );	
			}
			else
			{
				//for division
				verylong base = 0;
				verylong divisor = 0;

				//base 10 radix (decimal)
				zsread( "10", &base );

				//result = base^(diff-precision)
				if( (diff-precision) < 0 )
				{
					zsexp(base, 0 , &divisor);
				}
				else zsexp(base, (diff-precision) , &divisor);

				//divide to make numbers more manageable
				verylong trc_vl = 0;
				verylong crc_vl = 0;
				verylong quot	= 0;
				verylong rem	= 0;
				zdiv( total_result_count, divisor, &trc_vl, &rem);
				zdiv( current_result_count, divisor, &crc_vl, &rem);
		
				//if scaled 'total_result_count' (trc_vl) is still too big for double casting
				if( zcompare(trc_vl, max_long) != -1 ) 
				{
					boinc_fraction_done( 0.0 );	
				}
				//divide using scaled (by divisor) double vals
				else
				{
					boinc_fraction_done( ((double)ztoint(crc_vl)) / ((double)ztoint(trc_vl)) );
				}

				//free memory
				zfree( &trc_vl );
				zfree( &crc_vl );
				zfree( &quot );
				zfree( &rem );
				zfree( &base );
				zfree( &divisor );
			}			
		}

		//free memory
		zfree( &max_long );
}

#endif